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Abstract 

We present a mixed method for the linearized elas- 
ticity equations with independent approximation of 
the curl of the displacements. The curl can be seen 
as a drilling degree of freedom allowing for coupling 
with rotating objects and the direct application of 
moments of force. 

1 Introduction 

A drilling degree of freedom refers to a rotational 
degree of freedom for membrane elements, e.g. for 
linearized elasticity. Such degrees of freedoms refer 
to forces rather than displacements and are thus not 
natural to incorporate as enhancements of displace- 
ment fields. Nevertheless, the early attempts at for- 
mulating elements for planar elasticity with drilling 
degrees of rotations, these were seen as enhancements 
of the displacement field, see, e.g., Allman [31 [3] 
and Bergan and Felippa [3]. Such elements have to 
be carefully constructed using particular qualities of 
the basic elements they are meant to enhance, they 
arc difficult to generalize to higher order, and are 
plagued by stability problems. Since the drilling de- 
gree of freedom is a force-type variable, mixed meth- 
ods would seem more natural and were considered by 
Hughes and Brezzi [6], where several different meth- 
ods were proposed. The simplest of these was sub- 
sequently studied from a numerical point of view by 



Hughes, Masud, and Harari [7]. 

In this note, a different approach to mixed meth- 
ods for drilling degrees of freedom is taken for the 
case of isotropic linearized elasticity. Instead of ar- 
tificially adding the drilling degree of freedom to the 
energy functional, as in [7 , which can lead to numer- 
ical stiffening of the discrete problem, we introduce 
the drilling degree of freedom by splitting the weak 
form of the elasticity equations in a suitable way so 
that the curl of displacements is identified and can 
be replaced by an independent variable. 

The remainder is organized as follows: in Section 2 
we introduce the model problem and the split of the 
equations ; in Section 3 we remark on some differ- 
ent choices of finite element spaces; in Section 4 we 
compare and contrast our method with that of [7J; 
in Section 5 we present some numerical experiments, 
in particular comparing with [7J and with standard 
conforming methods. 



2 The Continuous Problem 

Consider a domain Q in R™ Bd , n s d = 2 or n s d = 3 with 
boundary T = Td U Tn, Td n Tn = 0, whose outward 
pointing normal is denoted n. The linear elasticity 



equations can be written 

-V-cr(w) = / infi 

u = onT D (1) 

n ■ cr(u) = g on Fn 

Here, with A and [i given material data, the stress 
tensor <x is defined by 



er(it) = 2it e(u) + A V • u 1 



(2) 



where it is the displacement field, 1 is the identity 
tensor, 

e(u) = - (V <g> u + (V ® w) T ) 

is the strain tensor, and / and g are given data. We 
have also used the notation 



(V-T) ( =£ 



J = l 



9a; i 



for the divergence of a tensor field t. Introducing the 
Hilbert space 

V = {v £ [H 1 ^)]* 1 ** :v = on T D }, 

the weak form of the elasticity equations is to find 
u E V such that 

/ 2/ie(u) : e(v) dx + / W-uV-vdx 
Jo. io 

= f ■ vdx + g ■ vds (3) 
Jn Jr N 

for all v € V. Here cr : e = a^Eij. 

The idea is to now rearrange the bilinear form in 
such a way that we can isolate the term 



/zV x it ■ V x v dx 



which corresponds to the rotational part of the form. 
It is also symmetric, leaving the remainder symmetric 
in turn. We shall then introduce a new variable p = 
V x it and construct a mixed method. Note that p 
is scalar in 2D. 



We remark that it is difficult to isolate the curl 
using the strong form of the equations, for instance 
using the identity V(V -it) = V- (V(g>it) T + V x Vxw, 
since the equations must be on the conservation form 
P]) in order to give the right weak boundary condi- 
tions in the corresponding weak form Q. Modifica- 
tions of the equations are thus best performed on the 
weak form after integration by parts. 

2.1 Rearranging the bilinear form 

In two dimensions we can explicitly write 
dui dv\ du2 dv2 



e(u) : s(v) 



dx\ dxi dx2 dx2 





( du 2 




f dv 2 


|- — ^ 




\dxi 


dx 2 J 


\dxi 


dx 2 J 



( du 2 


du! A 


( dv 2 


_dvi\ 


\dxi 


dx 2 j 


\dxi 


8X2 J 



and 

V x u ■ V x v = 
which means that 



/ \ / \ 1^ ^ du\ dvi 

e(u) : b(v) = -Vxm-Vxv + ■— — - 

2 axi ox\ 



du2 dv2 dui dv2 du2 dv\ 
8x2 dx2 0x2 dx\ dx\ 8x2 



or 



e(tt) : e{v) — -V xu-Vxo 



2^—i^-- ( 4 ) 



dxi dxi 



dx* dxi 



It is easily checked that relation (0| holds also in the 
three dimensional case. 

Introducing the variable p, we can write the 3D 
system on mixed form as the problem of finding 
(u,p) 6 V x Q, where 



Q = {v : v e [L 2 (n)} 3 }, 



such that 



a(u,v) + b(p,v) — f(v), Vv G V 



(5) 



2 



and 



where 



a(u, v) 



and 



b(q,u)-c(p,q)=g(q), Vg G Q (6) 



O I ZT^ ® X j ® X ' 



AV ■ uS7 ■ v dx, 



b(p,v) := / fip-Vxvdx, 



c{p, q) 



p ■ q dx, 



and by Helmholtz decomposition, cf . _5] , we split the 
force field / := V</> — V x h into a gradient field and 
a rotational field (with h x n = on Tn) and let 

f(v):= / V0-V(ia;+ / g-vds, 
Jn Jr N 



and 



.9(9) = / • <zdx- 



The split of / is not necessary but allows us to di- 
rectly apply distributed moments as external loads. 

By formally solving the second equation, we get 
p = (txV x u — h and from the first equation then 
formally follows 

-V • er(u) = V</> - V x h = f. 

In the 2D case the only difference is that p and h 
are replaced by scalar fields, Q = {v : «£ L2(£l)}- 



3 Finite element approxima- 
tion 

We replace the continuous spaces by discrete coun- 
terparts Q h C Q, V C V and pose the problem of 
finding (u h ,p h ) G V h x Q h such that 

a(u h ,v) + b(p h ,v) = /(»), Vvey' 1 (7) 



&(gy i )-c(p\g) = g(g), Vg g Q' 1 (8) 



The question then arises of what restriction there is 
on the combination of spaces V h and Q h with respect 
to stability of the discrete problem. 

A more common format for mixed methods is the 
case when c(-, •) = 0: find (u h ,p h ) G V h x Q h such 
that 



a(u h ,v) + b(p h ,v) = /(«), Vuey' 1 (9) 



and 



&(g,« h ) = s(g), Vg g Q h 



(10) 



This saddle point problem requires coercivity on the 
kernel of &(-, ■): there exists a constant a such that 



a(v,v) >a\\v\\ H i m \fv€K h , 



where 



K h = {v G F /l : 6(9, ») = Vg G Q fe }, 

and an inf-sup condition: there exists a constant /3 
such that 



inf sup — : 

q£Q h vev» \\v\\Hi(n)\\q\\L 2 (n) 



b(q,v) 



In such a mixed method, the space Q h cannot be cho- 
sen too large in order not to overconstrain the prob- 
lem, which would lead to stability problems. The 
presence of c(-, •) relaxes the inf-sup condition, but 
we still need coercivity on K h . Note that the bilin- 
ear form a(-, •) does not fulfill a Korn-type inequal- 
ity and is thus not coercive on the whole of V . In 
consequence, we need Q h to be large enough to get 
sufficient control of V x v to regain Korn's inequality 
on K h . In our case we will thus have the opposite 
problem of not letting Q h be too small in order to 
avoid stability problems. We note that the choice 

Q h := Vh x V h , 

i.e., functions in Q h are chosen as element- wise curl 
(Vh x ) of functions in V h , will give a method equiva- 
lent to the original elasticity problem, since p can 



3 



then be eliminated element- wise, and the mixed 
method will in that event be stable (for example a 
piecewise linear approximation in V and piecewise 
constant in Q ). More generally, if 

V h x V h C Q h 

the method will, for the same reason, be stable 
(though no gain comes from increasing the size of 
Q h ), but other choices will also work. For instance, 
in our numerical experience, the natural equal-order 
interpolation is stable. In the numerical examples we 
will show how some different choices of spaces behave. 

4 An alternative method 

The method of Hughes et al. [H [7] can, in the setting 
of using curl as an independent variable, be written in 
terms of minimization of a modified energy functional 

£ (u,p) = i / cr(u) : e(u) dfl 
2 Jn 

+ - / (p - V x u) ■ (p - V x u) dfl 
2 Jn 

- / f udn- / guds (11) 
Jn Jr N 

We note that here the new variable p is here intro- 
duced via a penalty-like functional with 7 as penalty 
parameter to be chosen. Setting 7 = /j, as recom- 
mended in [7] and minimizing the energy with re- 
spect to (u,p), followed by discretization, leads to 
the problem of finding (u h ,p h ) G V h x Q h such that 

/ (r(u h ) : e(v) dtt+ / x u h ■ V x vdVL 
Jn Jn 

-b(p h ,v) = /(«), VveV h (12) 

and 

b(q,u h )-c(p h ,q)=g(q), VqeQ h , (13) 

where we used the same split of the load vector as 
earlier. Here we can again choose Vh x V h C Q h and 
eliminate p h to obtain equivalence with the elasticity 
problem discretized in V h . Other choices have the 



drawback of leading to an additional stiffening of the 
discrete problem, as is obvious from the modified en- 
ergy functional The method is however coercive 
on the whole of V h so that in general we trade sta- 
bility for accuracy compared with ([7]) and ([S]) . In the 
following Section, we will compare the performance 
of the different methods for equal order interpolation 
on bilinear and linear elements. 

5 Numerical examples 

5.1 Convergence for different combi- 
nations of spaces 

We construct a right-hand side so that the exact so- 
lution is 

u {x 2 ~ l/2)x2Xl(l - Xl)(l - x 2 ) 

[ -X2Xl(xi — 1/2)(1 - Xl)(l - X2) 

We set A = /i = 1 and check the convergence in L 2 - 
norm of the Q1Q1 (equal order bilinear approxima- 
tions) and Q1P0 (bilinear displacements, piecewise 
constant curl) and compare with the standard bilin- 
ear method. In Fig. [1] we show the result which is 
second order convergence with slightly better error 
constants for the mixed methods. In Fig. [5] we show 
the convergence of the curl variable. Note that, at 
least on structured meshes, we get second order con- 
vergence of the curl for the Q1Q1 element. 

5.2 Comparison with the method of 
Hughes et al. 

We compare some different methods on an example 
consisting of a console defined by the domain < 
Xi < 1, < x 2 < 1) clamped at x\ = and with 
a surface traction g = (0,-1) at x 2 — 1; no volume 
load. In plane strain, with a Young's modulus of 
E = 1 and Poisson's ratio v = 0.3, the solution has 
the approximative "energy" 

\\u\\l~ [ tr(u) : e(u) dQ « 1.903697, 
Jn 

as given by Ainsworth et al. Q]. 
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In Fig. |3] we compare three different methods: 
mixed P1P1, standard PI (constant strain elements) 
and the method (fl"2")) - (fT3)) . which we call "Hughes 
method" in the following. We note that the new 
method is less stiff than the standard constant strain 
method, whereas Hughes method adds additional nu- 
merical stiffness. The same situation occurs if we 
take Ql — elements with equal order interpolation 
and with piecewise constant approximations for the 
rotation, see Fig. @] The equal order interpolated 
method is the least stiff. 

6 Concluding remarks 

We have introduced an approach to drilling degrees 
of freedom which is close to previously studied meth- 
ods [6l [7] but which introduces no artificial stiffening. 
The method, which is based on an independent ap- 
proximation of the curl of displacements, works best 
with equal order interpolation for displacements and 
curl. 
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Figure 2: Convergence of the curl. 
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Figure 3: Convergence of different linear methods in Figure 4: Convergence of different Qf-mcthods in 
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